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Abstract 

The  objective  of  this  article  is  to  demonstrate,  by  way  of  example(s),  how  to 
use  our  implementation  of  the  MLSL  method  for  model  independent 
parameter  estimation  to  calibrate  a  GSSHA  hydrologic  model.  The  purpose 
is  not  to  present  or  focus  on  the  theory  which  underlies  the  parameter 
estimation  method,  but  rather  to  carefully  describe  how  to  use  the  ERDC 
software  implementation  of  MLSL  that  accommodates  the  PEST  model 
independent  interface  to  calibrate  a  GSSHA  hydrologic  model.  Given  the 
computational  expense  associated  with  global  optimization,  we  will  initially 
consider  variations  of  our  MLSL  implementation  on  a  computationally 
efficient  test  problem  in  attempts  to  provide  the  interested  reader  with  an 
intuitive  sense  of  how  the  method  works. 


DISCLAIMER:  The  contents  of  this  report  are  not  to  be  used  for  advertising,  publication,  or  promotional  purposes. 
Citation  of  trade  names  does  not  constitute  an  official  endorsement  or  approval  of  the  use  of  such  commercial  products. 
All  product  names  and  trademarks  cited  are  the  property  of  their  respective  owners.  The  findings  of  this  report  are  not  to 
be  construed  as  an  official  Department  of  the  Army  position  unless  so  designated  by  other  authorized  documents. 

DESTROY  THIS  REPORT  WHEN  NO  LONGER  NEEDED.  DO  NOT  RETURN  IT  TO  THE  ORIGINATOR. 
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Preface 

This  report  describes,  by  way  of  multiple  examples,  how  to  use  model 
independent  software  to  calibrate  the  physics  based,  distributed 
parameter,  hydrologic  model  Gridded  Surface  Subsurface  Hydrologic 
Analysis  (GSSHA). 

Research  presented  in  this  technical  report  was  developed  under  the  U.S. 
Army  Corps  of  Engineers  Flood  and  Coastal  Storm  Damage  Reduction 
Research  and  Development  Program.  Dr.  William  Curtis,  Coastal  and 
Hydraulics  Laboratory  (CHL),  is  the  director. 

The  work  was  performed  by  Drs.  Brian  E.  Skahill  and  Charles  W.  Downer  of 
the  Hydrologic  Systems  Branch  (HF-H)  of  the  Flood  and  Storm  Protection 
Division  (HF),  U.S.  Army  Engineer  Research  and  Development  Center  - 
Coastal  and  Hydraulics  Laboratory,  and  Dr.  Jeffrey  S.  Baggett  of  the 
University  of  Wisconsin  -  La  Crosse.  At  the  time  of  publication,  Earl  V. 
Edris  was  Chief,  (HF-H);  Bruce  A.  Ebersole  was  Chief.  The  Deputy  Director 
of  CHL  was  Jose  E.  Sanchez  and  the  Director  was  Dr.  William  D.  Martin. 

COL  Kevin  J.  Wilson  was  the  Commander  and  Executive  Director  of 
ERDC,  and  Dr.  Jeffery  P.  Holland  was  the  Director. 
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1  Introduction 

Recent  research  at  the  U.S.  Army  Engineer  Research  and  Development 
Center  (ERDC)  has  focused  on  the  development  of  methodologies  or 
improvement  of  the  efficiency  of  native  algorithms,  for  the  computer-based 
calibration  of  hydrologic  and  environmental  models  (wherein  by  efficiency, 
we  mean  the  number  of  forward  model  calls  necessary  for  the  calibration 
algorithm  to  converge).  These  include,  among  others,  an  accelerated 
derivative-based  local  search  algorithm,  a  stochastic  global  optimization 
algorithm  for  intelligently  sifting  through  local  minima  to  find  a  global 
minimum  and,  most  recently,  a  state-of-the-art  evolutionary  strategy  for 
global  parameter  identification  of  difficult  problems  with  noise  or  other 
features  that  make  derivatives  estimation  difficult.  Minimizing  the  number 
of  required  model  runs  is  one  of  the  primary  factors  driving  the  research 
and  development  activities,  such  that  the  resulting  optimization  tool(s)  are 
more  compatible  with  the  computationally  expensive  physics-based  models 
that  are  becoming  more  commonly  used  within  the  practice  community. 

In  a  previous  technical  report  (Skahill  et  al.  2011)  we  recently  discussed 
how  to  use,  by  way  of  example,  our  implementation  of  the  Levenberg- 
Marquardt  (LM)  local  search  method  (Levenburg,  1944;  Marquardt,  1963), 
and  also  the  Secant  LM  (SLM)  method,  an  efficiency  enhancement  to  the 
LM  method,  for  computer-based  model  independent  hydrologic  model 
calibration.  The  context  for  this  article  will  be  directed  to  the  previously 
mentioned  stochastic  global  optimization  algorithm,  which  at  present  uses 
our  implementations  of  the  LM/SLM  methods  for  local  search.  The  LM 
method  has  features  that  make  it  attractive  for  model  calibration.  One 
feature  is  its  ability  to  readily  report  estimates  of  parameter  uncertainty, 
correlation,  and  (in)sensitivity  as  a  by-product  of  its  use  both  during  and 
after  the  parameter  estimation  process.  Another  feature  is  that  it  is  easily 
adapted  by  the  inclusion  of  various  regularization  devices  to  maintain 
numerical  stability  and  robustness  in  the  face  of  potential  numerical 
problems  (that  adversely  affect  all  parameter  estimation  methodologies) 
caused  by  parameter  insensitivity  and/or  parameter  correlation  (Menke, 
1984;  de  Groot-Hedlin  and  Constable,  1990;  Doherty  and  Skahill,  2006). 
Skahill  et  al.  (2009)  and  Skahill  and  Doherty  (2006)  both  provide  lengthy 
summaries  of  the  LM  method. 
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The  model  independent  LM  method  based  parameter  estimation  software 
PEST  (Doherty,  2004,  2007a,  2007b),  which  quantifies  model  to  measure¬ 
ment  misfit  in  the  weighted  least  squares  sense,  is  now  widely  used  to 
support  hydrologic  and  environmental  model  calibration.  In  addition  to  its 
traditional  groundwater  model  calibration  application  setting  (Zyvoloski  et 
al.,  2003;  Tonkin  and  Doherty,  2005;  Moore  and  Doherty,  2006;  Gallagher 
and  Doherty,  2007a),  it  is  now  employed  to  calibrate  ecological  models 
(Rose  et  al.,  2007),  land  surface  models  (Santanello  Jr.  et  al.,  2007)  and 
models  in  other  application  areas  including  nonpoint  source  pollution 
(Baginska  et  al.,  2003;  Haydon  and  Deletic,  2007),  surface  hydrology 
(Doherty  and  Johnston,  2003;  Gutierrez-Magness  and  McCuen,  2005; 
Kunstmann  et  al.,  2006;  Skahill  and  Doherty,  2006;  Doherty  and  Skahill, 
2006;  Gallagher  and  Doherty,  2007b;  Goegebeur  and  Pauwels,  2007;  Iskra 
and  Droste,  2007;  Kim  et  al.,  2007;  Maneta  et  al.,  2007),  and  surface  water 
quality  (Rode  et  al.,  2007). 

Skahill  et  al.  (2011)  focused  on  one  drawback  associated  with  LM-based 
model  independent  parameter  estimation  as  implemented  in  PEST;  viz., 
that  it  requires  estimates,  based  on  finite  differences,  of  the  derivatives  of 
the  objective  function  with  respect  to  the  model  parameters.  We  briefly 
discussed  the  secant  LM  method  and  then  presented  examples  which 
demonstrated  benefits  that  can  be  derived  in  terms  of  improved  inverse 
model  run-efficiency  when  using  the  SLM  method  rather  than  LM  for 
calibrating  a  GSSHA  hydrologic  model  deployment  for  the  Goodwin  Creek 
Experimental  Watershed  (GCEW). 

However,  as  mentioned,  another  drawback  of  the  LM  method  is  that  it  is  a 
local  search  method.  Thus,  if  there  are  different  “regions  of  attraction”  in 
parameter  space,  its  solution  will  lead  to  just  one  of  possibly  many  objective 
function  minima,  the  particular  one  that  is  found  is  dependent  upon  the 
user-supplied  set  of  initial  parameter  values.  Stochastic  global  optimization 
(GO)  can  be  employed  as  a  remedy.  Stochastic  global  optimization 
algorithms  estimate  the  global  minimum  of  the  objective  function  by 
initiating  local  searches  from  global,  randomly  sampled  points.  The  local 
and  global  phases  can  be  iterated  and/or  the  local  searches  may  be  initiated 
at  some  or  all  of  the  globally  sampled  points.  Stochastic  global  optimization 
algorithms  are  guaranteed  to  converge,  with  probability  one,  to  the  global 
minimum  as  the  sample  size  approaches  infinity.  Stronger  convergence 
properties  are  possible  for  some  stochastic  algorithms,  as  we  mention 
below.  Moreover,  probabilistic-based  stopping  criteria  can  be  developed  for 
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stochastic  global  optimization  methods  (Rinnooy  Kan  and  Timmer,  1987a, 
b;  Torn  and  Zilinskas,  1987);  however,  an  a  priori  computational  budget 
may  preclude  any  concern  regarding  termination  criteria. 

One  would  like  to  utilize  stochastic  GO  methods  that  are  not  only  reliable 
in  finding  the  global  minimum,  but  also  efficient  in  the  sense  that  they 
minimize  the  return  to  previously  visited  local  minima  in  parameter  space. 
A  modeler  would  possibly  also  like  to  receive  some  information  on  the 
locations  of  non-global  minima,  especially  if  these  minima  are  little 
different  in  magnitude  from  the  global  minimum,  but  are  widely  separate 
from  it  in  parameter  space.  With  the  understanding  that  efficient  and 
reliable  optimization  methods,  possibly  constrained  by  a  predetermined 
computational  budget,  that  are  capable  of  efficiently  finding  the  locations 
of  other  good  minima  in  addition  to  an  estimate  of  the  global  minimum, 
are  needed  to  identify  hydrologic  and/or  environmental  models,  Skahill  et 
al.  (2009)  implemented  a  more  efficient  and  reliable  stochastic  global 
optimization  algorithm  than  what  is  currently  available  in  PEST. 

Ideally,  we  would  like  to  perform  a  single  local  search  within  the  region  of 
attraction  of  each  local  minimum.  This  would  not  only  ensure  that  each 
local  minimum  is  identified  just  once,  but  also  that  in  fact  we  find  all  local 
minima.  But  we  also  want  to  employ  a  method  that  works  well  if  one  has  a 
predetermined  computational  budget  in  that  for  a  given  effort  it  compares 
favorably  with  other  methods.  Clustering  methods  were  designed  to 
accommodate  these  requirements.  They  are  variants  of  Multistart  (the 
Multistart  method  samples  points  from  a  uniform  distribution  over  the 
feasible  parameter  space  and  starts  a  local  search  from  each  of  the  sample 
points)  and  the  basic  idea  behind  them  is  to  group  close  points,  sampled 
from  the  feasible  parameter  space  and  for  which  the  specified  groups 
presumably  relate  to  actual  regions  of  attraction  in  parameter  space,  and  to 
apply  a  single  local  search  procedure  within  each  identified  cluster.  Either 
reduction;  wherein  sampled  points  associated  with  the  highest  objective 
function  values  are  temporarily  removed,  or  concentration;  wherein  the 
sampled  points  are  transformed  through  application  of  one  or  a  few 
iterations  of  a  local  search  procedure,  is  employed  to  identify  a  reduced 
sample  as  part  of  the  clustering  process  in  order  to  provide  some  assurance 
that  in  fact  the  specified  groups  correspond  to  regions  of  attraction  of  actual 
local  minima.  Clustering  methods  are  often  iterative  in  that  the  global  and 
local  phases  are  repeated  sequentially  until  a  stopping  criterion  is  satisfied. 
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With  clustering  methods,  it  is  possible  that  one  cluster  intersects  multiple 
regions  of  attraction;  hence,  the  global  minimum  could  be  missed,  or  that 
one  region  of  attraction  contains  more  than  one  cluster,  thus  allowing  for 
the  same  local  minimum  to  possibly  be  identified  more  than  once.  Multi 
Level  Single  Linkage  (MLSL)  is  a  clustering  method  that  was  developed  to 
reduce  the  probability  of  not  finding  a  local  minimum  or  of  finding  a  local 
minimum  more  than  once  (Rinnooy  Kan  and  Timmer,  1987a,  b). 

MLSL  mimics  clustering  by  calculating  a  critical  distance  rk  at  each 
iteration,  k.  This  critical  distance  can  be  used  to  build  clusters,  but  instead, 
in  MLSL,  the  decision  as  to  whether  a  local  search  is  to  be  initiated  from  a 
given  reduced  sample  point  is  simply  based  on  whether  there  exists 
another  reduced  sample  point  within  the  distance  rk  of  the  given  point 
with  a  corresponding  lower  objective  function  value.  The  critical  distance 
rk  is  reduced  at  each  iteration. 

Under  certain  assumptions,  MLSL  has  stronger  convergence  properties 
than  simpler  stochastic  global  optimization  algorithms.  First,  if  the 
algorithm  continues  forever,  the  number  of  local  searches  performed  is 
finite.  Second,  if  rk  tends  to  zero  with  increasing  k,  then  every  local 
minimum  will  be  identified  in  finite  time  with  probability  one. 

We  implemented  MLSL  and  it  uses  the  ERDC  implementation  of  the  LM  or 
SLM  method  for  local  search.  Our  MLSL  implementation  follows  that  of 
Rinnooy  Kan  and  Timmer  (1987a,  b)  with  a  slight  modification  to  some¬ 
times  avoid  repeatedly  finding  the  same  local  minima.  Our  MLSL 
implementation  is  summarized  in  Figure  1.  The  following  stopping  criteria 
were  included  as  part  of  our  MLSL  implementation: 

1.  The  objective  function  has  not  been  lowered  over  a  specified  number  of 
local  searches, 

2.  A  specified  maximum  number  of  local  searches  have  been  performed, 

3.  The  expected  number  of  minima,  in  the  Bayesian  sense,  exceeds  the 
number  of  identified  distinct  local  minima  by  less  than  0.5  (Rinnooy  Kan 
and  Timmer,  1987a,  b), 

4.  A  specified  maximum  number  of  MLSL  iterations  have  been  performed, 

5.  The  objective  function  has  not  been  lowered  over  a  specified  number  of 
MLSL  iterations. 
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Objective.  \  X  c  ft*. 

Initialization.  Let  S*  denote  the  set  of  identified  local  minima.  Start  with 
S*  =  0.  Let  the  integer  N  >  o,  y  €  (0,1]  and  <&,  <&  >  o  be  fixed,  let  k  :=  l  and 


Sampling  phase.  Use  a  uniform  distribution  to  sample  N  points  from  X 
at  step  k. 

Reduction  of  the  sample.  Sort  the  entire  sample  of  kN points  in  order 
of  increasing  objective  function  values.  Select  the  y  kN  points  with  the  low¬ 
est  objective  function  values.  This  resultant  set,  Rk,  is  called  the  reduced 
sample. 

Selection  of  the  set  Sk .  Apply  a  local  search  to  every  point  xi  in  the  re¬ 
duced  sample  set  which  satisfies  the  following: 

1.  There  is  no  other  point,  y,  in  Rk  or  S*  such  that  dix^y  and 

f(y)<f(xi ) 

2.  d[xi,S*)>d2 

3.  Alocal  search  has  not  already  been  applied  to  Xi. 

Stopping  phase.  If  the  termination  criteria  are  not  satisfied,  then  return 
to  the  sampling  phase  with  k  =  k  + 1. 

Figure  1.  Our  implementation  of  the  Multi  Level  Single  Linkage  (MLSL)  algorithm. 

An  instance  of  our  MLSL  implementation  stops  when  any  one  of  the  above 
criteria  are  satisfied.  It  can  be  utilized  by  simply  appending  the  MLSL 
search  parameters  to  the  end  of  the  control  data  section  of  a  working 
SLM_CHL  input  control  file  (see  Skahill  et  al.  2011).  Hence,  as  with  our 
independent  LM/SLM  implementations,  only  slight  modification  to  a 
working  PEST  input  control  file  is  required  to  also  use  our  MLSL 
implementation. 


Skahill  et  al.  (2009)  performed  an  efficiency  comparison  of  our  MLSL 
implementation  with  Trajectory  Repulsion  (Skahill  and  Doherty,  2006),  the 
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Shuffled  Complex  Evolution  (SCE)  general  purpose  global  optimization 
method  (Duan  et  al.  1992;  Duan  et  al.  1993)  and  the  Covariance  Matrix 
Adaption  Evolutionary  Strategy  (CMAES)  general  purpose  optimization 
method  (Hansen  and  Ostermeier,  2001;  Hansen  et  al.  2003),  using  an  eight 
parameter  Hydrological  Simulation  Program-FORTRAN  (HSPF)  (Bicknell 
et  al.,  2001)  hydrologic  model,  a  ten  parameter  Fast  All-season  Soil  Strength 
(FASST)  state-of-the-ground  model  (Frankenstein  and  Koenig,  2004),  and 
a  sixteen  parameter  Gridded  Surface  Subsurface  Hydrologic  Analysis 
(GSSHA)  (Downer  and  Ogden,  2003a,  2003b)  hydrologic  model.  In 
particular,  using  the  three  previously  mentioned  environmental  models, 
Skahill  et  al.  (2009)  compared  the  efficiencies,  in  terms  of  the  number  of 
model  calls  required  to  achieve  a  given  objective  function  value,  of  their 
implementations  of  Trajectory  Repulsion  and  MLSL  with  that  of  SCE  and 
CMAES  as  implemented  in  PEST  (Doherty,  2004,  2007a,  2007b).  To 
examine  their  relative  efficiencies,  thirty  trials  of  MLSL,  Trajectory 
Repulsion,  CMAES,  and  SCE  were  conducted  with  each  of  the  three  distinct 
environmental  model  structures.  While  PEST  includes  a  Trajectory 
Repulsion  implementation,  our  Trajectory  Repulsion  implementation  was 
utilized  so  that  the  local  search  implementation  (SLM)  would  be  identical  to 
that  employed  with  MLSL. 

Results  from  the  numerical  experiments  were  obtained  by  computing  the 
average  over  thirty  trials  of  the  best  objective  function  value  obtained  after  a 
specified  number  of  model  calls.  Skahill  et  al.  (2009)  pointed  out  that  they 
were  very  wary  of  comparing  different  software  packages  in  this  manner,  for 
a  package,  and  the  methodology  which  it  encapsulates,  always  performs 
best  when  operated  by  its  designers.  This  is  because  program  settings, 
particularly  those  pertaining  to  termination  and  convergence  criteria,  can 
have  a  huge  effect  on  the  performance  of  a  method;  a  non-expert  in  the  use 
of  a  particular  package  may  not  be  aware  of  the  optimal  settings  to  use, 
especially  in  difficult  cases.  Hence,  we  do  not  pretend  that  their  results 
provide  a  comprehensive  basis  for  assessment  of  the  comparative 
performance  of  SCE-PEST,  CMAES-PEST  and  MLSL.  We  hope,  however, 
that  they  do  provide  a  basis  for  at  least  a  “ball  park”  comparison  of  the 
methods  for  the  particular  calibration  cases  they  considered. 

Based  on  the  numerical  experiments  involving  thirty  trials  with  each  global 
optimization  method,  and  for  each  of  the  three  model  structures,  Skahill  et 
al.  (2009)  recommend  MLSL  over  Trajectory  Repulsion  for  environmental 
model  independent  LM -based  stochastic  global  optimization.  Moreover, 
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their  results  also  suggested  potential  utilization  of  MLSL  over  SCE  and 
CMAES,  except  for  the  case  where  only  a  very  limited  computational  budget 
is  available,  in  which  case  CMAES  might  be  preferable. 

The  objective  of  this  article  is  to  demonstrate,  by  way  of  example(s),  how 
to  use  our  implementation  of  the  MLSL  method  for  model  independent 
parameter  estimation  to  calibrate  a  GSSHA  hydrologic  model.  The 
purpose  is  not  to  present  or  focus  on  the  theory  which  underlies  the 
parameter  estimation  method,  but  rather  to  carefully  describe  how  to  use 
the  ERDC  software  implementation  of  MLSL  that  accommodates  the  PEST 
model  independent  interface  to  calibrate  a  GSSHA  hydrologic  model. 

Given  the  computational  expense  associated  with  global  optimization,  we 
will  initially  consider  variations  of  our  MLSL  implementation  on  a 
computationally  efficient  test  problem  in  attempts  to  provide  the 
interested  reader  with  an  intuitive  sense  of  how  the  method  works. 
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2  Examples 

The  steps  necessary  to  use  our  implementation  of  the  MLSL  GO  method 
will  now  be  demonstrated  and  documented  while  applying  it  to  initially 
calibrate  a  computationally  efficient  test  problem  and  subsequently  a 
GSSHA  hydrologic  model  for  the  Goodwin  Creek  Experimental  Watershed 
(GCEW)  (Senarath  et  al.  2000;  Downer  and  Ogden,  2003b).  The  general 
approach  simply  involves  slight  modification  of  an  existing  functional 
input  control  file  associated  with  the  ERDC  implementation  of  LM/SLM. 
The  interested  reader  is  referred  to  our  recent  report  for  details  on  how  to 
interface  a  given  forward  model  with  the  ERDC  LM/SLM  implementations 
(Skahill,  Downer,  and  Baggett,  2011). 

Example  1 

In  this  first  example,  we  will  document  the  steps  necessary  to  use  the  ERDC 
implementation  of  the  MLSL  stochastic  global  optimization  method  to 
calibrate  the  computationally  efficient  test  problem  defined  directly  below 
in  Figure  2. 


Prt*lil«*iii  I  n..!  ji\  <«ii'l  (<>  minimi/'': 

M \  .  *  .  *| 

^  f  1  MIH  ;i|  •  j*!  1  •  vin(  j»|  *  *«  D  •  -dtilrjlj 

ivlu'n  Ti  J.l  -  0.1  V  imil  £  0. 2, .[*2  G  [0.2;.  Flio  glulaiil  uiiu  :»  ft  I  Mir  viiv  by  «.yiuiui,tn'|. 

Thi*  |>n»bl«‘ui  Iuik  it  bin*  h  i*f  1»m.iI  minium.  l  id-  uwikr-  il  rlialJiMiuiut*  »•!-«'  im  im«t ImmI  will  ftml  rttlier  of  tin1  •?. b Ji.il 
inluiiu.i  ivi-ily 


Problem 


Figure  2.  Implementation  of  MLSL  with  many  local  minima. 
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Step  01  -  Obtain  forward  model 

The  C  source  code  provided  directly  below,  once  built  (it  was  named 
tp2.exe),  computes  model  generated  observations  for  a  given  parameter  set 
((pi,pa)).  The  program  reads  in  a  file  named  “tp2_parms.txt”,  a  file  with  two 
floating  point  data  values,  one  specified  on  the  first  row,  and  the  other 
specified  on  the  second  row  of  the  file,  respectively.  The  first  row  entry  is  the 
value  for  the  first  parameter,  p,;  whereas,  the  second  row  entry  is  the  value 
for  the  second  parameter,  p2.  The  87  model  computed  observations  are 
generated  and  written  to  a  file  named  “tp2_0ut.txt” ,  with  one  floating  point 
data  value  specified  per  line. 

////////////////////////////////////////////////////// 

//Example  Test  Problem  for  testing  MLSL  code 

// 

//Author:  Brian  E.  Skahill 

// 

///////////////////////////////////////////////////// 

#include  "stdio.h" 

#include  "conio.h" 

#  include  "stdlib.h" 

#include  "string.h" 

#include  "ctype.h" 

#  include  "math.h" 

#include  "time.h" 

void  main(  void ) 

{ 

inti; 

double  h[87],  value,  m[3],  x=o.o; 

FILE  Mata,  Mead; 

data  =  fopen(  "tp2_0ut.txt",  "w" ); 
read  =  fopen(  "tp2_parms.txt",  "r" ); 

////////////////////////////////////////////////////// 

//Read  the  parameter  values  -  m 

////////////////////////////////////////////////////// 

for(  i  =  1;  i  <=  2;  i++  ) 

{ 


} 


value  =  0.0; 

fscanf(  read,  "%lf\n",  &value  ); 
m[i]  =  value; 
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////////////////////////////////////////////////////// 

//Compute  h 

////////////////////////////////////////////////////// 

for(  i  =  o;  i  <=  86;  i++  ) 

{ 

x  =  3.1  +  o.i5*i; 

h[i]  =  (sin(m[i]*x)+sin(m[2]*x)); 

} 

////////////////////////////////////////////////////// 

//Write  the  computed  h  to  file 

////////////////////////////////////////////////////// 

for(  i  =  o;  i  <=  86;  i++  ) 

fprintff  data,  "%i6.i4E\n",  h[i] ); 

fclosef  data ); 

} 

Step  02  -  Collect  observed  data 

Observations  for  the  test  problem  can  be  generated  by  running  the  model 
executable  for  the  problem,  tp2.exe,  whose  source  code  was  presented  and 
briefly  described  in  the  previous  step,  with  input  parameter  values  set  to  p, 

=  1  andp2  =  2/3,  or  vice  versa.  Upon  execution  of  tp2.exe  at  pi  =  1  andp2  = 
2/3  (viz.,  the  global  minimum),  the  (synthetic)  observations  for  the  test 
problem  will  be  located  in  the  file  named  “tp2_0ut.txt” ,  as  mentioned 
above.  The  observations  for  the  test  problem  will  be  discussed  further  below 
and  clearly  noted  at  that  time. 

Step  03  -  Prepare  instruction  file  for  forward  model  output  file 

As  with  the  ERDC  LM/SLM  implementations,  our  MLSL  implementation 
is  model  independent.  Hence,  there  is  a  need  with  each  forward  model  call 
during  an  inverse  model  run  to  be  able  to  read  from  one  or  more  model 
generated  output  files  the  model  generated  observations  relevant  to 
objective  function  evaluation.  An  instruction  file  is  the  means  to  satisfy 
this  requirement.  As  was  also  mentioned  in  our  recent  technical  report, 
our  software  was  written  to  accommodate  the  popular  PEST  model 
independent  and  input  control  file  protocol;  hence,  the  interested  reader  is 
directed  to  Doherty  (2004,  2007a,  b,  c)  for  additional  details  regarding  the 
PEST  model  independent  interface,  and  in  particular,  instruction  file 
development.  The  model  instruction  file  for  the  test  problem,  named 
“tp2_out.ins” ,  is  provided  in  Appendix  1  (Please  note  that  all  of  the 
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appendices  associated  with  this  technical  report  are  available  in  a  separate 
document  made  publicly  available  on  the  GSSHA  Knowledge  Hub  at 
https://knowiedge.usace.armv.mii/l.  It  is  the  basis  for  reading  model  generated 
observations  from  the  model  ( tp2.exe )  output  file  named  “tp2_0ut.txt” . 

Step  04  -  Prepare  template  file  for  forward  model  input  file 

For  this  example,  the  model  input  file  named  “tp2_parms.txt”  contains  the 
two  input  parameter  values  that  are  designated  as  adjustable.  During  an 
inverse  model  run  (i.e.,  model  calibration),  this  model  input  file  will  need 
to  be  updated  every  time  there  is  a  need  to  evaluate  the  model  at  a  new 
location  in  adjustable  model  parameter  space.  To  support  the  interface  of 
the  model  for  this  test  problem  with  the  independent  PEST  and  ERDC 
LM/SLM  implementations,  a  PEST  template  file  was  prepared  for  this 
model  input  file.  The  contents  of  the  template  file  named  “tp2_parms.tpl” 
are  shown  directly  below. 

ptf$ 

$pi  $ 

$p2  $ 

Step  05  -  Prepare  preliminary  input  control  file 

We  are  now  ready  to  generate  an  input  control  file,  the  main  input  file  for 
execution  of  the  LM  method  not  only  associated  with  the  PEST,  but  also  our 
own  independent  LM/SLM  implementation.  The  path  that  we  take  for  this 
two  parameter  test  problem  will  be  different  from  the  path  that  we  pursued 
to  generate  an  input  control  file  for  the  GSSHA  GCEW  model,  as  we  docu¬ 
mented  in  our  recent  technical  report  (Skahill  et  al.  2011).  In  this  case,  we 
will  simply  take  an  existing  working  control  file  and  manually  modify  it  in  a 
text  editor  so  that  it  is  specific  to  our  two  parameter  test  problem  being 
considered  in  this  example.  In  particular,  we  will  manually  modify  the 
control  file  listed  in  Appendix  10  of  our  recent  technical  report  (Skahill  et  al. 
2011).  Relevant  changes  to  the  control  file  will  be  mentioned  so  that  the 
reader  can  interpret  and  follow  the  necessary  changes.  The  final  working 
input  control  file  for  this  example  problem,  named  tp2_ipst,  is  listed  in 
Appendix  2.  As  was  mentioned  in  our  recent  technical  report,  the  interested 
reader  is  referred  to  the  PEST  documentation  for  explanations  related  to 
control  files.  The  second  row  of  the  control  data  section  of  the  control  file  is 
changed  from  “16  233  16  o  1“  to  “2  8720  1”  because  we  now  have  two 
adjustable  model  parameters,  two  parameter  groups,  and  87  observations 
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rather  than  16  adjustable  model  parameters,  16  parameter  groups,  and 
233  observations.  Interpreting  the  parameter  groups  and  parameter  data 
sections  of  the  input  control  file  is  relatively  straight  forward,  and  the 
interested  reader  is  directed  to  the  PEST  documentation  and  also  our  recent 
report  (Skahill  et  al.  2011)  for  specific  details  related  thereof.  The  name  of 
the  one  observation  group  is  now  obsgroupi  rather  than  tmf.  The  87 
observed  data  points  are  listed  in  the  observation  data  section  of  the  input 
control  file,  and  they  are  all  equally  assigned  a  weight  of  one.  The  model 
command  line  section  of  the  control  file  lists  the  model  for  this  example 
and,  as  previously  noted,  it  is  tp2.exe.  And  the  two  relevant  input  and 
output  files,  also  previously  discussed,  are  listed  in  the  model  input/output 
section  at  the  very  end  of  the  control  file. 

Step  06  -  Verify  control  file  is  functional 

One  can  check  to  see  if  there  are  any  errors  with  the  input  control  file  as 
now  prepared  by  typing  the  following  at  the  command  prompt  and  pressing 
enter: 


PESTCHEK  tp2_i 

In  so  doing,  one  would  see  the  following  record  of  PESTCHEK  execution  at 
the  command  prompt  (Figure  3): 


C  \data\FVn_lilork\CHL-GSSHfl-MLSL\Example6\tp2>p@stchek  tp2_1 

PESTCHEK  Ueroion  12  1  Watermark  Numerical  Computing 

Errors  - > 

No  errors  encountered. 

Warnings  - > 

Consider  setting  JflCUPDATE  to  999  to  implement  Broyden's  Jacobian  updating. 
C  \data\FVH_ldorfc\CHL-GSSHA-MLSL\Exaople6\tp2> 


Figure  3.  Terminal  display  for  PESTCHEK  execution. 


Step  07  -  Modify  control  file  for  use  with  ERDC  MLSL  implementation 

To  employ  our  ERDC  MLSL  implementation,  two  additional  rows  of  input 
data  are  appended  to  the  end  of  the  control  data  section  of  the  control  file 
tp2_i.pst,  as  shown  below.  The  four  entries  on  the  first  row  appended  (now, 
the  next  to  last  row  of  the  control  data  section)  to  the  end  of  the  control  data 
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section  of  the  control  file  are  specific  to  our  independent  ERDC  SLM 
implementation,  and  the  description  and  function  for  each  of  the  four 
entries  is  provided  in  Step  to  of  Example  l  in  our  recent  report  (Skahill  et  al. 
2011).  So,  as  was  indicated  earlier  in  this  report,  with  our  MLSL 
implementation,  one  can  either  employ  LM  or  SLM  for  local  search.  The 
first  entry  on  the  second  row  appended  (now,  the  last  row  of  the  control 
data  section)  to  the  end  of  the  control  data  section  of  the  control  file  is  a 
floating  point  value  specifying  a  distance  threshold  that  is  used  for 
comparison  during  execution  of  our  implementation  of  MLSL.  In  particular, 
if  during  a  given  local  search  with  MLSL,  the  distance  computed  between 
the  current  location  in  adjustable  model  parameter  space  and  any  of  the 
previously  computed  parameter  upgrade  vectors,  obtained  either  during  the 
existing  or  with  any  of  the  previously  performed  local  searches,  is  less  than 
this  specified  threshold  value,  then  the  current  local  search  is  prematurely 
terminated  with  the  assumption  that  it  has  progressed  into  a  region  of 
attraction  of  a  previously  visited  local  minimum.  The  impact  of  this  input 
entry  has  not  been  explored  with  much  detail,  and  it  is  often  set  to  zero,  as  it 
is  specified  below.  Further  exploration  for  the  specification  of  this  input 
value  is  encouraged,  not  only  for  MLSL,  but  also  for  the  stochastic  GO 
method  Trajectory  Repulsion  (Skahill  and  Doherty,  2006).  The  subsequent 
four  entries  on  the  second  row  appended  (now,  the  last  row  of  the  control 
data  section)  to  the  end  of  the  control  data  section  of  the  control  file  are  all 
specific  to  our  ERDC  MLSL  implementation.  Relative  to  our  specification  of 
our  MLSL  implementation  outlined  in  Figure  1  above,  the  following  four 
input  values  are  N,  y,  o,  and  d2,  respectively.  N  is  an  integer  input  value; 
whereas,  the  remaining  three  are  floating  point  values.  N,  y,  and  o  impact 
the  efficiency  and  effectiveness  of  MLSL.  The  second  example  in  this  report 
will  explore  their  impact  using  the  computationally  efficient  test  problem 
from  this  first  example. 

*  control  data 
restart  estimation 
2  87  2  o  1 

I I  single  point  too 
5.0  2.0  0.3  0.03  10 
100000.0  3.0  0.001  o 
0.1  noaui 

30  0.00144  0.0014 

III 

o  o  1.0  o 
o  40  0.10  4  0.1 
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Step  08  -  Calibrate  test  problem  using  ERDC  MLSL  implementation 

Our  implementation  of  MLSL  was  employed  to  calibrate  the  computa¬ 
tionally  efficient  two  parameter  test  problem,  in  a  model  independent 
manner,  using  the  prepared  input  control  file  tp2_i.pst  (and  the  related 
input  files),  by  typing  the  following  at  the  command  prompt  and  pressing 
enter: 

chl_mlsl 

Upon  execution  of  our  current  implementation  of  MLSL,  the  user  is 
subsequently  prompted  for  additional  input  parameters.  Figure  4 
summarizes  its  complete  execution.  In  particular,  upon  execution  of  our 
implementation  of  MLSL,  the  user  is  asked  for,  in  sequence, 

1.  The  name  of  the  modified  input  control  file, 

2.  Whether  the  program  is  to  write  or  read  (thus  enabling  the  reproduction  of 
past  efforts)  the  file  named  “RandomNumberSeedsprn” ,  which  contains 
random  numbers  used  during  program  execution, 

3.  An  unsigned  long  that  provides  the  basis  for  computation  of  the  previously 
mentioned  random  numbers, 

4.  The  maximum  number  of  iterations  to  perform  during  MLSL  execution, 

5.  The  maximum  number  of  local  searches  to  perform  during  MLSL 
execution, 

6.  The  maximum  number  of  local  searches  to  perform  with  no  objective 
function  improvement, 

7.  The  maximum  number  of  MLSL  iterations  to  perform  with  no  objective 
function  improvement,  and 

8.  The  value  which  determines  whether  the  previous  and  current  computed 
best  objective  function  value  estimates  are  to  be  treated  as  effectively  the 
same,  or  distinct. 

The  MLSL  generated  output  file  that  will  be  of  most  interest  to  the  user  is 
named  “ct_mlsl.rec”,  and  its  contents  for  this  example  are  listed  in 
Appendix  3.  Examining  its  contents,  we  see  that  the  file  includes  an  echo 
of  the  MLSL  input  parameters,  and  then  for  each  MLSL  iteration, 

1.  A  listing  of  the  sampling  phase, 

2.  The  reduced  sample  set, 

3.  The  value  for  the  critical  distance  threshold  (named  alpha  in  the  file), 

4.  A  summary  of  the  treatment  of  each  member  of  the  reduced  sample  set, 
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and  then  the  file  ends  with  a  summary  of  the  stochastic  GO  run,  including 
the  number  of  local  searches  that  were  performed,  the  number  of  local 
searches  that  yielded  distinct  local  minima,  the  number  of  MLSL 
iterations,  the  minimum  objective  function  value  identified,  at  what  point 
during  MLSL  execution  the  minimum  objective  function  value  was 
identified,  the  total  number  of  forward  model  calls,  and  the  reason  for 
terminating  MLSL  execution.  For  example,  with  this  computationally 
efficient  two  parameter  test  problem,  MLSL,  with  the  input  as  specified 
above  and  also  listed  in  Appendix  3,  identified  the  global  minimum  in  the 
first  iteration  of  its  execution  with  the  first  local  search,  and  it  terminated 
execution  after  five  additional  MLSL  iterations  (and  one  additional  local 
search  that  yielded  another  local  minimum)  yielded  no  further  objective 
function  improvement.  One  hundred  and  seventy-six  forward  model  calls 
were  required  to  find  the  global  minimum. 


Figure  4.  Terminal  display  for  MLSL  method  (example  1). 


Example  2 

In  this  example,  we  will  further  examine  MLSL  execution  using  the 
computationally  efficient  two  parameter  test  problem  that  was  presented 
in  the  previous  example.  In  particular,  we  will  examine  the  impact  of  the 
input  parameters  N,  y,  and  o  on  the  performance  of  MLSL;  viz.,  its 
efficiency  (summarized  by  the  average  and  standard  deviation  of  total 
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model  calls  for  100  MLSL  trial  runs)  and  effectiveness  (summarized  by  the 
number  of  failures  to  find  the  global  minimum  among  too  MLSL  trial 
runs).  We  will  do  so  by  performing  too  trials  for  13  unique  N,  y,  and  o 
permutations.  Since  all  13  too  MLSL  trial  runs  are  very  similar  in  terms  of 
execution,  we  will  discuss  how  to  do,  by  way  of  example,  one  instance,  and 
then  summarize  the  results  for  all  13  too  MLSL  trial  runs  in  the  table  we 
present  below. 

Step  01  -  Prepare  100  MLSL  input  files  to  automate  execution  of  100 
MLSL  trial  runs 

To  enable  the  automatic  execution  of  too  MLSL  trial  runs,  too  MLSL 
input  files  were  prepared  whose  contents  are  like  what  is  shown  below. 

The  only  difference  between  each  separate  input  file  is  the  third  input 
value  within  each  file;  viz.,  the  unsigned  long  that  provides  the  basis  for 
computation  of  the  random  numbers  which  impact  the  sampling  phase 
which  occurs  in  each  iteration  of  MLSL  execution.  One  hundred  files  were 
prepared  with  contents  like  what  is  shown  below  and  the  too  files  were 
numbered  in  the  following  manner:  tp2_i_i.in,  tp2_i_2.in,  tp2_i_3.in, 

...,  tp2_i_99.in,  andtp2_i_ioo.in. 

tp2_l 

Y 

4063522982 

50 

too 

99 

5 

0.0025 

Step  02  -  Prepare  batch  file  for  execution  of  multiple  MLSL  runs 

A  batch  file  named  chl_mlsl_all.bat  was  prepared  and  its  contents  are 
listed  directly  below.  It  supports  the  execution  of  multiple  MLSL  runs. 

@ECHO  ON 

IF  [%i]==[]  GOTO  INPUTERROR 
Set  Count= 

:LOOP  (The  count  advances  one  each  time  through  the  LOOP} 
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::  Start  Counted  process. 

Set /A  Count +=  l 
@echo  %Count% 


chl_mlsl  tp2_i.pst  <  tp2_i_%Count%.in 

copy  RandomNumberSeeds.prn  RandomNumberSeeds_%Count%.prn 
copy  ct_mlsl.rec  ct_mlsl_%Count%.rec 
copy  tp2_i.rec  tp2_i_%Count%.rec 

IF  NOT  [%Count%]==[%i]  GOTO  LOOP 

Goto  END 

:INPUTERROR 

@ECHO  USAGE:  TPL  {Number  of  runs} 

@ECHO. 

:END 

Step  03  -  Calibrate  model  100  times  with  100  unique  calls  to  MLSL 

The  computationally  efficient  two  parameter  test  problem  was  calibrated 
too  times  by  way  of  too  MLSL  trial  runs  that  were  performed  by  typing 
the  following  at  the  command  prompt  and  pressing  enter: 

chl_mlsl_all  too 

Table  l  summarizes  the  results  from  all  13  separate  (calibration)  cases  of 
too  MLSL  trial  global  optimization  runs.  The  observed  differences  for 
each  case  are  completely  a  function  of  the  MLSL  input  parameters  N,  y, 
and  o. 

Table  1.  Results  from  all  13  separate  cases  of  100  MLSL  trial  global  optimization  runs. 
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30 
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10 
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10 

02 
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10 

01 
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20 
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0.05 

16 
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767  7 

775.93 
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706  47 
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92.64 
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100.62 
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103.65 
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Number  0#  fall  jr«  to  find  gloo*' 
minimum 

0 

4 

1 

1 

1 
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Example  3 

In  this  example,  we  will  describe  the  steps  necessary  to  perform  a  MLSL  GO 
run  to  calibrate  a  GSSHA  hydrologic  model  for  the  GCEW.  The  noted 
GSSHA  GCEW  hydrologic  model  is  the  same  model  that  was  also 
considered  in  the  examples  in  our  recent  report  which  documented  how  to 
use  our  independent  implementations  of  the  LM/SLM  local  search  methods 
for  model  calibration  (Skahill  et  al.  2011).  To  that  end,  before  continuing 
any  further  with  this  example,  wherein  we  will  discuss  how  to  use  the 
stochastic  GO  method  MLSL  to  calibrate  the  GSSHA  GCEW  hydrologic 
model,  the  active  reader  is  first  directed  to  complete  the  first  and  third 
examples  in  our  recent  report  (Skahill  et  al.  2011).  In  particular,  the  point  of 
departure,  relative  to  completion  of  example  1  and  example  3  from  our 
recent  report,  and  our  intent  to  use  MLSL,  is  to  simply  augment  the  control 
data  section  of  the  input  control  file  from  example  3  in  our  recent  report 
(Skahill  et  al.  2011);  viz.,  the  input  control  file  gc_i_bu2.pst,  as  shown 
below.  The  four  new  additions,  which  are  relevant  to  MLSL  execution,  are 
highlighted  yellow,  and  the  one  edit,  which  most  likely  will  effectively  not 
allow  there  to  be  a  full  update  of  the  model  Jacobian  during  any  of  the  local 
searches,  is  highlighted  orange. 


*  control  data 
restart  estimation 
16  233  16  o  1 
1 1  single  point  too 
5.0  2.0  0.3  0.03  10 
5.05.0  i.oe-3 
0.1  noaui 


30  .00544  .0054 
111 

1  o  100000000000.0 


o  20  0.05  4  0.1 


o 


After  the  slight  alteration  to  the  input  control  file  noted  above,  all  that  is 
now  required  is  to  call  MLSL  at  the  command  prompt 


chi  mlsl 


The  inputs  that  were  further  specified  at  the  command  prompt  for  MLSL 
execution  are  shown  below  (Figure  5).  And  the  contents  of  the  file 
ct_mlsl.rec,  which  summarize  the  salient  aspects  of  the  MLSL  GO  run  to 
calibrate  the  GSSHA  GCEW  hydrologic  model,  are  listed  in  Appendix  4. 
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Figure  5.  Terminal  display  for  MLSL  method  (example  3). 


Example  4 

A  second  calibration  run  for  the  GSSHA  hydrologic  model  for  the  GCEW 
was  also  performed  using  our  MLSL  implementation,  this  time  with  the 
MLSL  input  parameter  o  set  equal  to  one  rather  than  four,  as  in  the 
previous  (third)  example,  in  attempts  to  require  more  exploration  of  adjust¬ 
able  model  parameter  space.  The  termination  criteria  for  this  example, 
relative  to  the  previous  (third)  example,  were  specified  at  the  command 
prompt  in  a  manner  to  encourage  a  more  exhaustive  search  than  what  was 
performed  in  the  third  example  (Figure  6).  The  contents  of  the  file 
ct_mlsl.rec,  which  summarize  the  salient  aspects  of  the  MLSL  GO  run  to 
calibrate  the  GSSHA  GCEW  hydrologic  model,  are  listed  in  Appendix  5. 
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Figure  6.  Terminal  display  for  MLSL  method  (example  4). 
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3  Results  and  Discussion 

Example  2  explored  the  impact  of  the  MLSL  input  parameters  N,  y,  and  o 
on  overall  algorithm  efficiency  and  effectiveness  by  considering  thirteen 
separate  cases  of  100  MLSL  global  optimization  trial  runs  which  were 
directed  at  calibrating  the  computationally  efficient  two  parameter  test 
problem  that  was  originally  presented  in  example  1.  Examining  the  table 
that  was  presented  in  example  2  (Table  1),  which  summarized  the  thirteen 
separate  numerical  experiments  involving  100  trials  in  each  case,  the 
following  observations  are  given: 

1.  With  y  and  o  fixed,  as  N  increases,  so  does  the  average  for  the  total  number 
of  forward  model  calls. 

2.  With  N  and  y  fixed,  as  o  decreases,  we  observe  that  the  average  and 
standard  deviation  for  the  total  number  of  forward  model  calls  both 
increase,  modestly,  and  also  that  the  effectiveness  of  the  MLSL  algorithm 
improves.  For  example,  we  see  that  with  N,  y,  and  o  equal  to  20,  0.05,  and 
16,  respectively,  we  have  2  failures  for  the  100  trials,  but  only  1  failure  for 
the  100  trials  when  o  =  8,  and  no  failures  for  the  100  trials  when  o  =  4,  2, 
or  1.  The  MLSL  algorithm  is  more  effective  with  the  lower  o  values  because 
the  critical  distance  which  is  used  to  characterize  the  size  of  the  clusters  is  a 
function  of  o  (see  Figure  1).  At  a  lower  o  value,  cluster  size  is  smaller; 
hence,  the  opportunity  for  a  single  cluster  to  contain  multiple  regions  of 
attraction  is  reduced. 

3.  With  N  and  o  fixed,  as  y  increases  so  does  the  average  and  standard 
deviation  for  the  total  number  of  forward  model  calls,  simply  because  y 
controls  the  potential  number  of  local  searches  to  be  performed  during 
each  iteration  of  MLSL  execution,  and  as  y  increases  so  does  that  potential. 

Based  on  our  own  experience  to  date,  we  recommend  that  users  of  our 
implementation  of  MLSL  start  with  a  value  of  one  for  o,  and  choose  values 
for  N  and  y  that  take  into  consideration  the  expense  associated  with  a 
single  forward  model  call  and  the  understanding  that  MLSL  is  an  iterative 
stochastic  GO  method.  And  we  recommend  that  termination  criteria  be 
specified  in  a  manner  consistent  with  the  level  of  importance  assigned  to 
finding  the  global  minimizer  for  the  given  application. 
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In  examples  3  and  4  we  used  our  implementation  of  MLSL  to  calibrate  the 
GSSHA  hydrologic  model  for  the  GCEW  which  was  originally  calibrated  as 
part  of  the  work  effort  documented  in  our  recent  report  using  our  indepen¬ 
dent  LM/SLM  local  search  method  implementations.  For  examples  3  and 
4,  in  consideration  of  the  expense  associated  with  a  single  forward  model 
call  for  the  GSSHA  GCEW  hydrologic  model,  we  used  the  SLM  method  for 
local  search,  in  such  a  manner  that  no  measures  were  employed  to 
mitigate  against  the  potential  that  our  approximation  to  the  true  model 
Jacobian  may  become  poor  after  some  optimization  iterations.  In  example 
3,  517  total  model  calls  were  required  to  calibrate  the  model  using  MLSL 
with  the  given  input  parameters  (as  specified  above  and  also  at  the  very 
beginning  of  Appendix  4).  And  in  example  4,  with  the  value  for  o  reduced 
relative  to  its  value  used  for  example  3  (from  4  in  example  3  to  one  in 
example  4),  and  also  using  stricter  termination  criteria,  3832  total  model 
calls  were  required  to  calibrate  the  GSSHA  GCEW  model  using  MLSL  with 
the  given  input  parameters  (as  specified  above  and  also  at  the  very 
beginning  of  Appendix  5).  As  a  result,  in  example  3,  7  local  searches  were 
performed  and  7  distinct  local  minima  were  identified;  whereas,  in 
example  5,  55  local  searches  were  performed  and  55  distinct  local  minima 
were  identified.  With  example  3,  the  minimum  and  maximum  identified 
objective  function  values  were  58.126056  and  70.299196,  respectively,  and 
the  minimum  value  was  found  from  the  fourth  of  the  seven  total  local 
searches.  With  example  4,  the  minimum  and  maximum  identified  object- 
tive  function  values  were  57.825337  and  77.795517,  respectively,  and  the 
minimum  value  was  found  from  the  fiftieth  of  the  fifty-five  total  local 
searches.  Figures  7  and  8  are  plots  of  the  transformed  observed  and 
simulated  flows  that  constitute  the  minimum  objective  function  values 
that  were  identified  in  Examples  3  and  4,  respectively.  The  figures  provide 
the  interested  reader  with  a  means  to  effectively  compare  the  minimum 
objective  function  values  that  were  obtained  in  the  two  examples  in  terms 
of  a  statistical  summary  (i.e.,  R2  and  Nash-Sutcliffe  efficiency  score)  and 
also  by  way  of  a  scatter  plot  of  the  transformed  observed  data  and  its 
model  simulated  counterparts.  Clearly,  there  is  little  difference  in  terms  of 
fit  with  the  observed  data  between  the  models  associated  with  the  mini¬ 
mum  objective  function  values  that  were  identified  in  Example  3  and  4. 
The  salient  difference  is  that  the  minimum  identified  in  Example  4 
required  approximately  an  order  of  magnitude  more  model  calls.  Possibly 
even  of  greater  interest  to  the  engaged  reader  would  be  to  compare  the 
results  presented  in  Figures  7  and  8  with  the  results  that  were  displayed  in 
a  similar  manner  for  Example  3  in  our  recent  report  (Skahill  et  al.  2011), 
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wherein  the  same  GSSHA  GCEW  hydrologic  model  was  calibrated  using 
the  independent  ERDC  SLM  implementation,  employed  in  the  same 
manner  as  it  was  for  the  local  searches  within  MLSL  in  Examples  3  and  4. 
Performing  such  a  comparison,  one  would  see  comparable  results  were 
obtained  when  using  SLM  simply  once,  with  no  cyclic  updating,  in  terms 
of  the  fit  achieved  between  the  observed  data  and  their  model  simulated 
counterparts,  wherein  an  associated  final  objective  function  value  of  59.56 
was  obtained  in  only  62  total  model  calls. 


Conparison  of  Heasured  Data  oith  Sinulated  Data 


Figure  7.  Plot  of  transformed  observed  and  simulated  flows,  associated  with  the  minimum 
objective  function  value  identified  in  Example  3. 
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Conparison  of  Heasured  Data  oith  Sinulated  Data 


Figure  8.  Plot  of  transformed  observed  and  simulated  flows,  associated  with  the  minimum 
objective  function  value  identified  in  Example  4. 
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4  Summary 

In  this  report,  we  demonstrated  how  to  use  our  implementation  of  the 
stochastic  global  optimization  technique  called  Multi  Level  Single  Linkage 
(MLSL),  which  uses  our  LM/SLM  method  implementations  for  local  search. 
We  considered  four  separate  examples  to  demonstrate  its  practical  use.  All 
four  examples  emphasized  how  simple  it  is  to  use  our  MLSL  implementa¬ 
tion  after  a  given  forward  model  has  already  been  interfaced  with  our 
LM/SLM  software  for  local  search  (see  our  recent  report  (Skahill  et  al. 
2011)).  The  first  example  discussed  all  of  the  steps  necessary  to  use  our 
MLSL  implementation  to  calibrate  a  computationally  efficient  two 
parameter  test  problem.  The  second  example  also  involved  calibrating  the 
computationally  efficient  two  parameter  test  problem  that  was  considered 
in  example  1.  It  considered  a  series  of  13  separate  numerical  experiments, 
each  involving  too  unique  MLSL  runs  for  a  given  set  of  input  parameters,  in 
attempts  to  provide  the  active  reader  with  an  intuitive  feel  for  how  the 
principal  MLSL  input  parameters  affect  MLSL  algorithm  efficiency  and 
effectiveness.  The  third  and  fourth  examples  both  involved  demonstrating 
how  to  calibrate,  using  our  MLSL  implementation,  the  same  GSSHA  GCEW 
hydrologic  model  that  we  also  considered  in  our  recent  report  (Skahill  et  al. 
2011)  wherein  we  demonstrated  howto  use  our  LM/SLM  local  search 
implementations  for  GSSHA  hydrologic  model  calibration.  The  fourth 
experiment  was  designed,  among  others,  to  provide  an  opportunity  to 
examine  MLSL  termination  criteria  and  their  impact  on  overall  algorithm 
execution.  We  provided  some  initial  guidance  regarding  the  specification  of 
MLSL  input  parameter  values  and  termination  criteria. 

Our  quick  comparison  between  the  results  that  were  obtained  in  Examples 
3  and  4,  and  also  with  the  SLM  local  search  calibration  run  that  was 
performed  and  documented  upon  in  our  recent  report  (Skahill  et  al.  2011) 
encouraged  us,  based  on  our  limited  applications  to  date,  to  recommend 
that  GSSHA  hydrologic  modelers  currently  use  our  SLM  implementation, 
possibly  with  prior  information,  rather  than  MLSL  for  GSSHA  hydrologic 
model  calibration.  However,  if  needed  and/or  desired,  then  we 
recommend  MLSL  as  the  global  optimization  tool  to  use  to  calibrate 
GSSHA,  rather  than  the  SCE  method,  which  is  currently  employed  for 
automatic  calibration  of  GSSHA  models  (please  see  Skahill  et  al.  (2009) 
for  further  details).  The  previously  mentioned  comparisons  underscore  the 
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need  to  provide  for  an  effective,  efficient,  and  adaptive  means  for 
calibrating  spatial  hydrology  models  such  as  GSSHA.  We  see  that  to  be 
possible  by  merging  ideas  from  previous  contributions,  in  particular,  by 
merging  ideas  from  Skahill  et  al.  (2009)  and  Doherty  and  Skahill  (2006). 
This  is  our  current  research  and  development  focus  for  GSSHA  hydrologic 
model  calibration. 

The  user  of  our  ERDC  implementation  of  the  MLSL  method  accepts  and 
uses  the  software  at  his/her  own  risk.  Any  questions,  comments,  and/or 
concerns  regarding  the  use  of  our  ERDC  software  implementation  of  the 
MLSL  method  with  the  GSSHA  model  should  be  directed  to  the  first 
author.  The  interested  reader  is  encouraged  to  contact  the  second  author 
with  any  questions,  comments,  and/or  concerns  related  to  the  GSSHA 
hydrologic  simulation  model. 
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